### PPP
### Descriptive tables and figures

rm(list = ls())

library(foreign)
library(xtable)
library(Hmisc)
options(scipen=999)

setwd("~/Dropbox/Pollution and Public Perceptions in China/Data and Analysis/Data")

# loading the data
load("ppp_cleaned")

# change wd for output (tables)
setwd("~/Dropbox/Pollution and Public Perceptions in China/Manuscript/Tables")

p <- p2

# creating vector of covariates
colnames(p)

p$educo <- as.numeric(as.character(p$educo))
p$nghbr_trust <- as.numeric(as.character(p$nghbr_trust))
p$kids <- as.numeric(as.character(p$kids))
p$c_sat <- as.numeric(as.character(p$c_sat))
p$l_corr <- as.numeric(as.character(p$l_corr))



vars <- c("gender","ageo","hukou","educo","incomeo","soe","ccp","insider","kids","married")

# recoding SOE
p$soe[p$soe==2] <- 0
mean(na.omit(p$soe))

# then create dfs that are the only the selected vars
stats <- p[vars]

# for table, dividing income by 10k - calling it incomeo
mean(na.omit(stats$income))

stats$Income <- (stats$income)/10000

# getting rid of old income var and replacing w/ new (in same order)
stats$incomeo <- NULL
stats <- stats[c(1,2,3,4,10,5,6,7,8,9)]
colnames(stats)

colnames(stats) <-  c("Female", "Age", "Urban Hukou (0 No / 1 Yes)", "Education (1-8 level)", "Annual Income(10,000 RMB)", "SOE Employment (0 No / 1 Yes)", "CCP Member (0 No / 1 Yes) ", "Regime Insider (0 No /1 Yes)", "Kids(0 No / 1 Yes)","Married (0 No / 1 Yes)")

Obs. <- apply(stats, 2, function(x){sum(!is.na(x))})
Means <- sapply(stats, mean, na.rm=TRUE)
Median <- sapply(stats, median, na.rm=TRUE)
SD <- sapply(stats, sd, na.rm=TRUE)
Min <- sapply(stats, min, na.rm=TRUE)
Max <- sapply(stats, max, na.rm=TRUE)

# binding these columns
summ <- round(cbind(Obs.,Means,Median,SD, Min, Max),digits=2)
# inverting hukou value to avoid confusion
View(summ)
summ <- xtable(summ,digits=c(0,0,2,0,2,0,0))

## TABLE A2
####
sink("~/Dropbox/Pollution and Public Perceptions in China/Manuscript/Tables/summ.tex")
print(summ,floating=FALSE,font.size="tiny",latex.environments="center")
sink()

#####
## wd for figures 
setwd("~/Dropbox/Pollution and Public Perceptions in China/Manuscript/Figures")

### vertical
## FIGURE A3
pdf(file = "bars_vertical.pdf", height = 8, width = 6)
par(mfrow =c(3,2))
barplot(prop.table(table(factor(p$l_sat,levels=1:5))), xlab="Local Govt. Satisfaction",ylab="Proportion",main="Local Gov. Satisfaction")
barplot(prop.table(table(factor(p$c_sat,levels=1:5))), xlab="Central Govt. Satisfaction",ylab="Proportion",main="Central Gov. Satisfaction")
barplot(prop.table(table(factor(p$l_watch,levels=1:5))), xlab="Demand for Oversight of Local Gov.",ylab="Proportion",main="Demand for Local Oversight")
barplot(prop.table(table(factor(p$c_watch,levels=1:5))), xlab="Demand for Oversight of Central Gov.",ylab="Proportion",main="Demand for Central Oversight")
# barplot(prop.table(table(factor(p$econ,levels=1:5))), xlab="Confidence in the Chinese Econonomy",ylab="Proportion",main="Confidence in the Economy")
barplot(prop.table(table(factor(p$anti_west,levels=1:5))), xlab="Agreement w/ Anti-West Pollution Propaganda",ylab="Proportion",main="Agreement w/ Ant-West Prop.")
dev.off()

### horizontal
pdf(file = "bars_horizontal.pdf", height = 5, width = 9)
par(mfrow =c(2,3))
barplot(prop.table(table(factor(p$l_sat,levels=1:5))), xlab="Local Govt. Satisfaction",ylab="Proportion",main="Local Gov. Satisfaction")
barplot(prop.table(table(factor(p$c_sat,levels=1:5))), xlab="Central Govt. Satisfaction",ylab="Proportion",main="Central Gov. Satisfaction")
barplot(prop.table(table(factor(p$l_watch,levels=1:5))), xlab="Demand for Oversight of Local Gov.",ylab="Proportion",main="Demand for Local Oversight")
barplot(prop.table(table(factor(p$c_watch,levels=1:5))), xlab="Demand for Oversight of Central Gov.",ylab="Proportion",main="Demand for Central Oversight")
# barplot(prop.table(table(factor(p$econ,levels=1:5))), xlab="Confidence in the Chinese Econonomy",ylab="Proportion",main="Confidence in the Economy")
barplot(prop.table(table(factor(p$anti_west,levels=1:5))), xlab="Agreement w/ Anti-West Pollution Propaganda",ylab="Proportion",main="Agreement w/ Ant-West Prop.")
dev.off()

#######
#######
## Perceptions Table
######
######

#create unhealthy dummy
p$high <- 0
p$high[p$aqio>150]<- 1 

# no covariates
fit1 <- lm(day_pollute ~ high + parade + national, data = p)
m.vcov1 <- cluster.vcov(fit1, p$day)
se1 <- coeftest(fit1, m.vcov1)

# covars
fit2 <- lm(day_pollute ~ high + parade + national + hukou + insider + soe + kids + married + gender + ccp + ageo + income, data = p)
m.vcov2 <- cluster.vcov(fit2, p$day)
se2 <- coeftest(fit2, m.vcov2)

fit.list <- list(fit1,fit2)
ses <- list(se1[,2],se2[,2])
covars <- c("Unhealthy Pollution (Dummy)","Parade Period","Holiday","Hukou","Regime Insider","SOE Employee","Children","Married","Female","CCP","Age","Income")

sink("~/Dropbox/Pollution and Public Perceptions in China/Manuscript/Tables/high.perception.tex")
stargazer(fit.list, font.size = "large", digits = 4, 
          column.labels = "Perception that Today's Pollution is Not a Serious Problem", 
          covariate.labels = covars,
          dep.var.labels.include = FALSE,
          model.numbers = FALSE,
          style = "qje", no.space = FALSE, keep.stat = c("n"),
          star.cutoffs = c(.10, .05, .01),
          se = ses,
          header = F, notes = "robust standard errors clustered by day in parentheses", float = F)
sink()

#### Table showing AQI(Continuous) as Cause of Day Pollute

# no covariates
fit3 <- lm(day_pollute ~ aqi + parade + national, data = p)
m.vcov3 <- cluster.vcov(fit3, p$day)
se3 <- coeftest(fit3, m.vcov3)

# covars
fit4 <- lm(day_pollute ~ aqi + parade + national + hukou + insider + soe + kids + married + gender + ccp + ageo + income, data = p)
m.vcov4 <- cluster.vcov(fit4, p$day)
se4 <- coeftest(fit4, m.vcov4)

fit.list <- list(fit3,fit4)
ses <- list(se3[,2],se4[,2])
covars <- c("AQI","Parade Period","Holiday","Hukou","Regime Insider","SOE Employee","Children","Married","Female","CCP","Age","Income")

# Table A4
sink("~/Dropbox/Pollution and Public Perceptions in China/Manuscript/Tables/aqi.perception.tex")
stargazer(fit.list, font.size = "large", digits = 4, 
          column.labels = "Perception that Today's Pollution is Not a Serious Problem", 
          covariate.labels = covars,
          dep.var.labels.include = FALSE,
          model.numbers = FALSE,
          style = "qje", no.space = FALSE, keep.stat = c("n"),
          star.cutoffs = c(.10, .05, .01),
          se = ses,
          header = F, notes = "robust standard errors clustered by day in parentheses", float = F)
sink()